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Abstract 

We report a study of exchange interactions in bulk Cr02 calculated from first principles. We 

a ■ 

g ■ considered three near neighbor Cr-Cr exchange interactions: the interaction between corner and 

\ body center atoms mediated through a single oxygen atom; the interaction between a Cr and the 

Cr directly "above" it in the (001) direction, also mediated by a single O atom; and the interaction 

o : 

between a Cr and its neighbor in the (100) direction, mediated by two intervening oxygen atoms. 
The interactions were calculated by rotating the moments of one or more of the Cr ions while 
j — . ' constraining the others to remain parallel. We then fit the resulting energy vs. angle data to the 

Heisenberg model and extracted exchange energy parameters with a least-squares method. We also 
calculated the exchange interactions using a "spin-spiral" technique, in which a relative angular 

ON . 

displacement was imposed upon Cr moments in adjacent cells. Similar results were obtained 



X 



with both approaches. The calculated T = K exchange interactions were subsequently used to 
determine the magnetization as a function of temperature via low-T spin-wave dispersion and a 
Monte-Carlo method. Reasonable agreement with experiment was obtained. 
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I. INTRODUCTION 



&O2 is one of only a few known ferromagnetic oxides, and is the only material which has 
been experimentally shown to be a ferromagnetic "half- metal"- 1 - 2 . A half- metal is a material 
which is a metal for one spin channel and an insulator for the other. &O2 crystallizes in the 
rutile crystal structure (Figured]), as do TiC>2, VO2, M11O2, RuC>2, and SnC>2. The existence 
of isostructural oxides with a variety of different electronic and magnetic properties makes 
the rutile system interesting for theoretical investigations of spintronics because one can 
envisage the growth of layered devices with the same crystal structure throughout. Since 
Cr0 2 offers such special opportunities for understanding oxide spintronics, it is important to 
establish how well our standard electronic structure tools work in dealing with the electronic 
and magnetic structure of this material. It is well known that they encounter difficulties 
in dealing with many transition metal oxides, including the very similar oxide VO2, which 
DFT also predicts to be a half-metal at OK,- but is observed to be an insulator. An addi- 
tional motiviation for understanding exchange interactions in Cr02 is the fact that its Curie 
temperature (T c = 386.5 K)&& is sufficiently close to room temperature that its magnetic 
properties are significantly degraded at room temperature, hindering potential spintronics 
applications. A better understanding may point the way to improvement. 



II. ELECTRONIC STRUCTURE OF Cr0 2 WITHIN DENSITY FUNCTIONAL 
THEORY 

The electronic structure and density of states were calculated using density functional 
theory^ (DFT). Our calculated density of states is similar to previous calculations. 1 The 
electronic structure of Cr02 can be understood by comparing the DOS of Cr02 with that 
of Ti02- Rutile Ti02 is a non-magnetic insulator. In the rutile structure, as in many 
transition metal oxides, each transition metal atom is surrounded by six oxygen nearest- 
neighbors. Each oxygen atom has three transition metal (TM) nearest neighbors. The TM 
ions form a body-centered tetragonal lattice, with a - ratio of approximately |. The oxygen 
ions form rows collinear with the Cr ions. The rows containing the body-center Cr ions can 
be considered to run in the (110) direction, in which case the rows containing the corner Cr 
ions run in the (110) direction. In addition to the lattice constants (a = b and c), a single 
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internal parameter is sufficient to describe the structure. This parameter may be taken to 
be the distance in the a-b plane between a Cr ion and its neighboring oxygen ions. The 
oxygen atoms in the top and bottom faces of the tetragonal cell form a rectangle with a 
Cr ion at the center. The - ratio and the internal parameter usually conspire such that 
the environment of each Cr ion is an approximate octahedron of oxygen ions. However, the 
octahedron is distorted because the square base of the octahedron with the Cr ion at its 
center is actualy a rectangle. In addition, the distance from a Cr ion to the oxygen ions at 
the apices of the octahedron may differ from the distance to the four equatorial oxygen ions 
that form the rectangular base. As shown in Figured! the edges of distorted octahedra form 
"ribbons" along the (001) direction. 

It is straightforward to show that if we treat this system in a tight-binding approximation 
in which the TM atoms only interact directly with the oxygen atoms (i. e. hopping matrix 
elements only connect nearest neighbors), there will be an energy gap separating the oxygen 
p-states and the TM d-states. The gap extends from the O-p onsite energy to the O-d onsite 
energy. This gap is apparent in Ti02, for which the oxygen p-states are filled and the Ti 
d-states are empty (Fig. [2]). When an energy gap occurs at the Fermi energy, it contributes 
significantly to reducing the energy of the structure, because all occupied states are pushed 
down in energy, while all unoccupied states are pushed up. In Cr0 2 , there are two additional 
electrons per TM atom compared to Ti02, so some of the rf-states above the gap must be 
occupied. 




FIG. 1: Rutile structure. For CrC>2, we use a = 4.42 and | 0.670 (experimental parameters). 

Let us consider three possibilities for the occupation of the d-states in Cr02- The first 
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FIG. 2: Density of States for rutile Ti02 

possibility is that the system remains non- magnetic, similar to Ti02 (Fig. [3]). In this case 
each Cr atom will have one spin-up and one spin-down <i-electron. From Fig. [3] we can see 
that there remains a gap between the O-p states and the Cr-d states, but it is significantly 
smaller than for TiC>2. If the interactions are strictly nearest neighbor connecting only O 
ions with Ti or Cr ions, then the O-p on-site energy in a tight-binding model would be at 
the top of the O-p states and the transiton metal d on-site energy would be at the bottom of 
the d-states. It should be understood that both states below the gap and the states above 
the gap are hybrids of O-p and TM-d. However, the lower set of states are predominantly 
O-p, and there are 3 of them per spin channel per O ion, and the higher lying states are 
predominantly TM-d, and there are 5 of them per spin channel per TM ion. 
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FIG. 3: Density of States for Nonmagnetic Cr02- 
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A second possibility is that the system becomes ferromagnetic, the experimentally- 
observed ground state. In this case there will be 1 + x majority spin electrons and 1 — x 
minority spin electrons (Fig. H]). In this spin configuration, energy is gained from the intra- 
atomic "Hund's first rule" energy, which in most versions of DFT is primarily a potential 
proportional to the local magnetization. Additional energy can be saved if x = 1, because 
this restores the gap for the minority spin channel. It is clear from Fig. H] that the system 
chooses x = 1 for the ferromagnetic case. The system does this by lowering the on-site en- 
ergy of the majority of Cr <i-states and raising the on-site energy of the minority Cr cf-states. 
If we were confident that there were only nearest neighbor interactions, the majority density 
of states would imply that the Cr majority d on-site energy was essentially degenerate with 
the O-p on-site energy. In reality, there are next-nearest-neighbor interactions that spread 
the O-p states above the O-p on-site energy and spread the Cr-d states below the Cv-d onsite 
energy. 
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FIG. 4: Density of States for Ferromagnetic Cr02- 

A third possibility is to maintain two electrons of the same spin on each Cr atom, but 
to make the system antiferromagnetic — for example, by having opposing moments on the 
corner and body-centered TM atoms of the 6-atom primitive cell (Fig. [5]). In this case, the 
spin-dependent DOS of the corner and body-center Cr ions are mirror images of each other. 
Assuming that the corner Cr ion moment is "up," its d on-site energy for spin up is pushed 
down while its d on-site energy for spin down is pushed up. The opposite will happen for 
the body-center Cr ion. The result is that there is no gap in the total DOS for either spin 
channel, but if one looks only at the DOS for one of the Cr ions, one can see an approximate 
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gap in one of the spin channels. Thus for this anti-ferromagnetic configuration, there would 
be a relatively large density of states at the Fermi energy for a given spin channel, but only 
on alernate Cr ions. 




FIG. 5: Density of States for Antiferromagnetic Cr02- 



Comparing these three possibilities, it is not surprising that the ferromagnetic state has 
the lowest energy, the nonmagnetic state the highest (1.02 eV above ferromagnetic) with 
the anti-ferromagnetic intermediate between the two (0.30 eV above ferromagnetic). Thus, 
the tendency to form a moment in Cr02 is very strong, and the energy associated with 
the ferromagnetic alignment of moments based on this initial test is moderately large within 
DFT. It should be recognized that other more complicated spin arrangements (e. g. , different 
antiferromagnetic states) may have lower energy than the simple one calculated here. 



III. EXCHANGE INTERACTIONS IN Cr0 2 

In order to investigate interatomic exchange interactions in Cr02 in more detail, we 
have calculated the near-neighbor exchange interactions along the (100), (001), and (111) 
directions by rotating moments within specially-constructed supercells. We fit the resulting 
relationship between the energy of the system and the angle of rotation to the Heisenberg 
model 

hi 

where |/x| = g^sS = 2/ig is the spin moment, g is the electron spin ^-factor, S — 1 is the 
spin number, and hb is the Bohr magneton. To make contact with the standard Heisenberg 

6 



model, we can pull the magnitude of the spin moment (2/ie) into the value of J and treat 
the spins as unit vectors. 

In addition to this supercell approach, we have taken advantage of a recently developed 
feature in the Vienna Ab-initio Simulation Package 6 (VASP) to calculate a so-called heli- 
magnetic state in which the moment in the n th magnetic layer is canted by an angle n6 with 
respect to the th layer. In so doing, we are able to calculate several orders of J n of the form 



via Fourier analysis. 

A. Near Neighbor Exchange Using Supercells 

All of the calculations in this study were performed within DFT 7 in the generalized 
gradient approximation- (GGA) and in the local (spin) density approximation with onsite 
Coulomb interactions (LSDA+U)^ using the Dudarev method,— for which we use U — J = 
3 eV. We perform all calculations using the VASP software 6 - and pseudopotentials gener- 
ated by Kresse et al.— . To calculate the near- neighbor exchange interactions, we created a 
supercell containing two rutile unit cells (using both experimental and DFT-relaxed lattice 
parameters), stacked in either the (100) (Fig. [6]) or (001) (Fig. [TJ) direction as appropriate. 
In all of the following calculations, we use an energy cut-off of 500 eV. For cells stacked 
along the (100) direction, we use a 5 x 9 x 15 Monkhorst-Pack 12 grid of k-points, a 9 x 9 x 7 
grid for supercells stacked along (001), and a 9 x 9 x 15 grid for the 6-atom cell used in the 
spin-spiral calculations. We also make use of the spin interpoloation method of Vosko-Wilk- 
Nusair.— Each of the 12-atom supercells has four Cr ions, whose magnetic moments we can 
individually constrain within the calculation. We chose three distinct magnetic configura- 
tions, designed to probe the exchange coefficients. In the first configuration, we rotated the 
moment of a corner Cr atom and held all other moments fixed using the constraining field 
method in VASP. In the second, we rotated the two Cr moments in the centers of their 
respective unit cells, and in the third we rotated a corner atom and its nearest center atom. 
A summary of the configurations can be found in Table [B 

To ensure that the systems would be sufficiently well-behaved, we rotated the moments 
through small angles (up to 60°). We fit the energy vs. angle data to A(l—cos6)+B, where A 




(2) 



n 
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FIG. 6: The (100) supercell, with Cr ions numbered for comparison to Tabled The bottom half 
of the structure in Figured] has been projected onto the (010) plane. 




FIG. 7: The (001) supercell, with Cr ions numbered for reference. The left half of the structure in 
Figure [Q has been projected onto the (010) plane. 



is equivalent to the combined exchange parameter J and B is simply the angle-independent 
component of the energy. 

For a given choice of supercell orientation, we have the following system of equations for 
the supercell method: 

^Case 1 — 8Jm + 2J 10 o/001 (3) 
^Case 2 = 16 J m (4) 
^Case 3 = 8 Jm + 4Jioo/001 (5) 

Using a least-squares technique for overdetermined systems of equations,— we can write 



AJ = b 



(6) 



C ri Cr 2 Cr 3 Cr 4 



Case 1 fixed rotated fixed fixed 



Case 2 rotated rotated fixed fixed 



Case 3 rotated fixed rotated fixed 



TABLE I: Magnetic configurations used to calculated exchange coupling. 



A T AJ = A T b 



(7) 




(8) 



a= A-J-b 



(9) 



where J is the calculated J column vector, a is the error in the fit, and 
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We summarize the calculations performed within GGA and LSDA+U for experimental 
and relaxed lattice parameters using both the supercell method in Table [Til Throughout this 
work, the terms "experimental" and "relaxed" as they appear in this table denote structures 
with the experimental and the GGA- or LSDA+U-relaxed lattice parameters, respectively. 

The results of the calculations for the three cases are summarized as follows: in each case, 
we find a near-perfect fit to the cosine function, provided that we restrict the fit to small 
angles (less than or equal to 60°), as we did with the original calculations. We can see clearly 
in Table [TfJ the anisotropic nature of the exchange, which is to be expected given the shape 
of the cell. Most interestingly, we find that the interaction between Cr neighbors along the 
(100) or (010) directions (parallel to the a or b axes) is antiferromagnetic. However, the 
strength and multiplicity of the other interactions is sufficient to lead to a ferromagnetic 
ground state. Considering the dependence on lattice parameter, we notice that the (001) 
and (111) interactions seem to be almost unchanged with the small (0.6%) change in lattice 
constant. Somewhat surprisingly, however, the (100) interaction (calculated within the 
GGA) increases (becomes more positive) by more than an meV under this small expansion 
of the lattice. We also note that the LSDA+U calculations predict a positive, though nearly 
vanishing, Ji o- 
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GGA 



LSDA+U {U 



J = 3 eV) 



z Experimental Relaxed Experimental 



Relaxed 



Jioo (meV) 4 -11.8 ± 2.5 -10.4 ±0.7 1.4 ±0.5 



1.1 ± 1.0 



Jooi (meV) 2 33.8 ± 5.6 33.8 ± 5.0 32.8 ± 0.3 



31.5 ±0.3 



J m (meV) 8 23.2 ±6.1 22.9 ± 5.0 24.5 ± 0.4 



24.6 ± 1.0 



TABLE II: Summary of all calculated exchange energies obtained using the supercell method. 
Uncertainties given arise from the error in the least-squares fit. In Jm, there is some (usually 
negligible) contribution to the error from the standard deviation of the values obtained through 
(100)- and (OOl)-stacked supercells. The quantity z is the coordination number for the given 
interaction — the number of such interactions acting on a given Cr ion. Note that the (100) and 
(010) directions are equivalent and are referred to as (100) throughout this work. The columns 
"experimental" and "relaxed" refer to the experimental or DFT-relaxed structure. 

B. Helimagnetism 

Helimagnetism is a noncollinear magnetic state in which the spins in adjacent layers along 
a certain direction are rotated with respect to one another by a fixed angle. Rutile MnC>2, for 
example, has been shown to exhibit helimagnetic ordering in the ground state.— We do not 
suspect that Cr02 is a helimagnetic material, but by setting up a helimagnetic spin state, 
we can investigate the exchange using a different approach. The recently-added spin spiral 
capabilities of VASP 14 allow us to calculate arbitrarily long-range exchange interactions 
within bulk CrC^. 

The spin spiral method modifies the periodic boundary conditions of the supercell ap- 
proach, imposing helimagnetic order on the magnetic structure as determined by the prop- 
agation vector q. The vector q and the angle between any two spins are given by 



where the azimuthal angle 6 is restricted to f. We choose the unit vector to be either 
the (100) or (001) direction, and allow £ to vary between and 1. Clearly, when £ = 0, we 
recover the ferromagnetic state. 



= q • *j 



(11) 




(12) 
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Because the unit cell contains two magnetic ions, varying the angles between neighboring 
layers of CrC>2 requires that one modify both £ and the initial orientation of the magnetic 
moments. For example, in this work, we use £ = 1 in conjunction with an antiferromagnetic 
configuration of the two moments in the unit cell to create a periodic magnetic system 
wherein the moment of each atom's nearest neighbor is antiparallel it. To obtain a system 
in which neighboring magnetic "layers" are oriented at an angle of | from one another, we 
use £ = \i so that each cell after the initial one is rotated by |. We then set up the moments 
in the initial cell such that the corner and body-centered Cr moments are oriented at the 
desired angle of |, leading to a smooth spin wave in the desired direction. 

In this work, we choose a relatively short, commensurate spin wavelength in order to 
simplify the analysis, although the method allows for more general configurations as well. 
Using different values of q, and thus different values of 8, we create a q spectrum. We then 
use Fourier analysis to extract the J n . These J n differ in meaning from the Js calculated 
using the supercell method; for example, J\ = 8Jm and J2 = 2Ji o/ooi- 

To calculate the helimagnetic state, we used a supercell composed of a single rutile unit 
cell. The angle of each subsequent Cr ion with respect to the first is given by ( fl2l) . After ac- 
quiring five points (including the zero-frequency point q=0) of the E(q) curve, we performed 
a discrete Fourier transform to obtain the first 4 J n . We find good agreement between the 
J\ calculated with (100) and (001) spin spirals, as expected. We also find a difference in 
sign between J2 in the (100) and (001) cases, in agreement with the larger supercell calcu- 
lations. Moreover, this method yields the additional parameters J3 and J4, corresponding 
to J211/112 and J200/002, respectively. These higher-order energies are smaller than the first- 
and second-order exchange energies, and will be neglected in further analysis. The results 
of the calculations are summarized in Table IIHI 

IV. COMPARISON WITH EXPERIMENTS 
A. Spin Wave Stiffness 

In an effort to benchmark our calculations against known experimental results, we 
have calculated the spin wave stiffness constant for Cr02 using the expression derived by 
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GGA LSDA+U 
Experimental Relaxed Experimental Relaxed 



•^100 


(meV) 


-12.0 


-12.2 


-1.3 


-1.8 


-Tool 


(meV) 


27.5 


29.8 


29.6 


26.1 


Jin 


(meV) 


20.8 


20.7 


25.6 


25.0 



TABLE III: Summary of calculated exchange interactions (in meV) using the spin spiral method 
(compare with Table [TTJ) . The effect of the change in lattice parameter is decidedly smaller (perhaps 
negligible in most cases) in the spin spiral method. Note that in this method the LSDA+U 
calculation also predicts a negative Jiooi albeit an order of magnitude smaller than in the GGA. 

Schlottmann:— 

Dioo = 2(J m + Ji 00 )Sa 2 (13) 

Axn = 2(J m + J oi)Sc 2 (14) 

where a and c are the lattice spacings in the appropriate directions and S is the spin number 
(1 for CrC^). These expressions can be easily understood as anisotropic extensions of results 
obtained for magnons in a one-dimensional chain (for which D = 2JSa 2 ). In his work, 
Schlottmann considers the spins as quantum operators such that Sj • Sj = S(S + 1), and 
he keeps this value separate from J. Additionally, he neglects Jioo in his expression for 
.Dioo- However, we use classical spins of magnitude 2/i# (although the units are collapsed 
into the exchange constant J as previously explained). Consequently, we must scale our Js 
by 1/ |a*| 2 = 1/4 in order to apply this expression. Further, our calculations indicate that 
Jioo it is roughly of the same order as J m and J oi, so we have included it in our analysis. 
Using this model, we calculate -Dioo an d -D oi for the various cells, exchange-correlation 
approximations, and methods considered throughout this work. Table HVl reviews the values 
we obtained. Experimentally, we find several mostly-consistent values for the spin wave 
stiffness obtained through different methods. All of the experimental values assume an 
isotropic stiffness constant. Ji et al. 20 fit the M(T) curve in order to obtain the coefficient 
on the T 3 / 2 term, from which they determine D = 1.8 x 10~ 40 Jm 2 . Zou et al— used 
magnetic force microscopy to determine the length and width of domain walls in Cr02, from 
which they were able to calculate D = 2.62 x 10~ 40 Jm 2 . Further, Rameev et alM used 
ferromagnetic resonance to measure the bulk magnon modes and obtained Db = 3 x 1CU 10 
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Oe cm 2 , which is equivalent to D = 0.57 x 10 40 Jm 2 via the relation D B = 2A/n M s ,— 
which is a bit smaller than other reported values. 

GGA LSDA+U 
Experimental Relaxed Experimental Relaxed 



Supercell .Dioo 


(10 
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io- 


-40 


J 


m 2 ) 


1.81 


1.96 


4.06 


3.94 


-Dooi 


(10 


X 


10" 


-40 


J 


m 2 ) 


3.91 


3.87 


3.91 


3.79 


■D aV g 


(10 


X 


10" 


-40 


J 


m 2 ) 


2.34 


2.46 


4.00 


3.89 


Spin Spiral .Dioo 


(10 


X 


10" 


-40 


J 


m 2 ) 


1.38 


1.35 


3.80 


3.56 




(10 


X 


io- 


-40 


J 


m 2 ) 


3.29 


3.46 


3.76 


3.45 


D aV g 


(10 


X 


10- 


-40 


J 


m 2 ) 


1.84 


1.84 


3.79 


3.52 



TABLE IV: Comparison of the calculated spin stiffness constants D for different methods of first- 
principles calculation. Here, D avg = (-DiooVAxn) 2 ^ 3 - 

Using our calculated Ds, we can predict the low-temperature spin-wave contribution to 
the magnetization as a function of temperature. A relatively straight-forward generalization 
of the argument found in Kitted for a cubic system allows one to write the spin-wave 
dispersion relation for small excitations and long wavelengths as 

u(k,k z ) =Dk 2 + D z kl k 2 = k 2 x + k 2 , D = D 100 , D z = D m (15) 

Integrating over a surface of constant u in fc-space, one obtains a density of states given by 

= ^d7W^ (16) 



Using this expression and the Planck distribution, we can calculate the coefficient B in the 
T 3 / 2 model 

M(T) = M(0)(1 - BT 3/2 ) (17) 
= 0.0587 1 1 fc3 /2 = 0.0587 V fe3/2 , , 

SQ 2S(J W0 + Jm) y/2S( Jooi + Jm) B SQ D\fD~ z B 1 } 

where Q is the number of magnetic ions per unit cell (2, in this case), V is the volume of 
the cell, and ks is Boltzmann's constant. Fitting the experimental 20 M(T) curve yields B = 
5 x 10~ 5 K~ 3//2 . Using the spin-wave stiffnesses shown in Table HVl we have, for supercells, 
B e G f A = 2.40 x 10- 5 K- 3 / 2 , B^ GA = 2.27 x 10~ 5 K" 3 / 2 , B e L f DA+u = 1.07 x lO" 5 K" 3 / 2 , and 
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b lsda + u = im x 10-5 K ~ 3/2 - For the s P in s P iral approach, B%g A = 3.43 x 1CT 5 FT 3 / 2 , 
B r G % A = 3.48 xlO" 5 K" 3 / 2 , B e L f DA+u = 1.16 xHT 5 K^ 2 , tmdB r L % A+u = 1.27xlCr 5 K^ 2 . 
Thus, the coefficient obtained from GGA is within a factor of two , while that derived from 
LSDA+U is off by about a factor of four. Assuming that DFT overestimates each exchange 
energy equally, this implies that our calculated values of J may differ from experimental 
values by about 50% for GGA and a factor of about 2.5 for LSDA+U (with U — J = 3 eV). 
In each case, the spin spiral numbers are slightly closer to experiment. Figure M shows the 
low-T M(T) curves from the calculated spin- wave dispersion compared to that from a fit to 
the experimental M(T) curve. 
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FIG. 8: The low-temperature M(T) curve. The GGA Spin Spiral and LSDA+U Supercell curves 
represent the extremes of the range of calculated M(T) curves. We compare against the actual 
experimental data£ and a low-T fit to these data. 



B. Curie Temperature 

In light of the favorable agreement between calculated and experimental spin stiffness, we 
subsequently attempted to calculate the magnetic ordering temperature of Cr02, comparing 
a mean field prediction to Monte Carlo simulations. A mean-field model using the calculated 
exchange parameters yields a Curie temperature several times larger than the measured value 
of 386.5 K.-£ The mean-field expression is given by 

k B T = ~J tot (19) 
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where J tot is equivalent to half of the energy difference between a ferro- and an antiferromag- 
netic configuration in a 6-atom (2-Cr) cell. Using this expression, we obtain a mean-field 
Curie temperature for Cr02 of 1160 K or 1240 K for the experimental and DFT-relaxed 
lattice parameters in the supercell method, respectively. This is somewhat surprising given 
the above analysis of our estimation of the exchange. However, it is not sufficient to consider 
only the low-temperature behavior, especially when the system in question has a magnetic 
ordering temperature above room temperature. In order to gain a simple yet ideally illumi- 
nating picture of the temperature dependence, we utilized a Monte Carlo simulation using 
the Metropolis-Hastings algorithm 24 with random numbers generated using the Mersenne 
Twister method 2 - 5 . For this simulation, we used a cubic grid of 10 x 10 x 10 unit cells, where 
a unit cell consists of a corner and body-centered Cr ion. Only Cr ions are considered, and 
they are treated as simple constant-magnitude magnetic moments that initially point along 
the z axis. Our calculations indicate that the constant-magnitude approximation should be 
valid as long as the angle between adjacent moments is less than 100°. 

We begin with a random spin configuration with the spin vectors chosen to be uniformly 
distributed on the unit sphere. In the Metropolis method, an iteration consists of a randomly 
chosen Cr ion being assigned a magnetic moment in a random direction. This will result in a 
change in energy AE from the old configuration. If AE is negative, meaning the new energy 
is lower, the new direction for that moment is kept. Otherwise, the new direction still has a 
probability of e _AjB//fcsT of being kept in its new orientation to simulate thermal agitation. If 
neither condition for keeping the moment's new direction is met, then the change is undone, 
and the lattice of spins remains unmodified until the next iteration. 

The calculation of AE at each step considers all nearest neighbors along (100), (010), 
(001), and (111) directions, using a Heisenberg interaction between moments with the cal- 
culated exchange constants for GGA and LSDA+U with experimental and DFT-relaxed 
lattice parameters. Figure M shows the simulated results for the magnitude of the net 
magnetization versus temperature compared to reported values.- When interpreting these 
data, one must must be cognizant of the fact that the Monte Carlo simulations exhibit 
several shortcomings — namely, that it will necessarily not be able to predict the correct 
low-temperature T-dependence (as it uses a classical model), that there exists an unphysical 
tail on the curve arising from finite-size effects in the lattice, and that we assume that ex- 
change remains constant with temperature, likely leading to an overestimation of the Curie 
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temperature. The errors in the shape of the curve at low temperature should not have an 
impact in the accuracy of the result, as each value of ksT is run independently. Further, 
the high-temperature tail can be accounted for assuming we know where to look for it. The 
remaining discrepancy, that the exchange will reduce in strength as temperature rises, is 
a limitation of exploring this behavior with first-principles calculations. Considering these 
issues, we can attempt to scale the experimental curve to match the calculated one, as we do 
in Fig. [91 By matching the portions of the curves at lower temperatures, we can determine 
with some confidence where the curve ought to go to zero and thus the true prediction of the 
model. By so doing, we find that the Monte Carlo simulation predicts a Curie temperature 
of between 450 K (for the GGA spin spiral method) and 650 K (for the LSDA+U supercell 
method), which, even in the latter case, is in much better agreement with experiment than 
the mean-field model. 
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FIG. 9: The calculated M(T) behavior predicts a Curie temperature of 450-550 K, depending on 
which exchange parameters one chooses from among those presented in this work. We compare 
with the experimental data 4 . 



V. CONCLUSIONS 

We have calculated the near neighbor exchange interactions for bulk Cr02 in the (100), 
(001), and (111) directions, finding them to be —12.0 meV, 24.5 meV, and 23.6 meV re- 
spectively. From our calculated spin stiffness parameters and the results of our classical 
Heisenberg Metropolis method, we obtain some confidence that DFT and VASP can deter- 
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mine the exchange interaction of &O2 with reasonable accuracy. Examining the calculated 
exchange parameters, we find that the sign of J100, both in the supercell and the equiva- 
lent spin spiral calculations, indicates the possibility of non-collinear behavior in CrC>2 if 
the exchange parameters are modified. Thus, a mixed interface between &O2 and another 
material (such as RUO2) might lead to non-collinear spins if the ratio between nearest and 
next-nearest neighbor interactions is pushed into a "favorable" zone. We investigate this 
possibility explicitly for Cr02-Ru02 interfaces in a companion paper by K. Chetry et alr§-. 
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